#include "stream.h"

static FILE *fh_load, *fh_save;
static char runname[80], load_name[80], save_name[80];
static char param[160], value[160], line[160];
static fftw_complex *tmp, *tmp2;
static double dx, dy, *x, *y, kx, ky;

void load_params(char** input, ctrl_ptr ptr)
{
if (input[1] == NULL)
  {
    printf("Usage:\t%s input\n",input[0]);
    exit(1);
  }
  strcpy(runname, input[1]);
  strcpy(ptr->runname, runname);
  sprintf(load_name,"%s.in",runname);
  sprintf(save_name,"%s.out",runname);
  fh_load = fopen(load_name,"r");
  if (fh_load == NULL) 
  {
    printf("%s\tno such file or file corrupt.\n",load_name);
    exit(1);
  }
  while (fgets(line, 160, fh_load) != NULL)
  {
    sscanf(line, "%s\t%s", param, value);
      if (strcmp(param,"#nx") == 0) ptr->nx = atoi(value);
      if (strcmp(param,"#t_skip") == 0) ptr->tskip = atof(value);
      if (strcmp(param,"#ny") == 0) ptr->ny = atoi(value);
      if (strcmp(param,"#lx") == 0) ptr->lx = atof(value);
      if (strcmp(param,"#ly") == 0) ptr->ly = atof(value);
      if (strcmp(param,"#max_time") == 0)   ptr->t_max = atof(value);
      if (strcmp(param,"#cfl_num") == 0)    ptr->cfl = atof(value);
      if (strcmp(param,"#U") == 0)    ptr->U = atof(value);
      if (strcmp(param,"#a") == 0)    ptr->a = atof(value);
      if (strcmp(param,"#ic_name") == 0)    strcpy(ptr->ic_name, value);
      if (strcmp(param,"#forcing") == 0)    ptr->FFlag = atof(value);
      if (strcmp(param,"#av_power") == 0)   ptr->p = atoi(value);
      if (strcmp(param,"#alpha") == 0)   ptr->alpha = atof(value);
      if (strcmp(param,"#a_force") == 0)   ptr->a_force = atof(value);
      if (strcmp(param,"#k_force") == 0)   ptr->k_force = atof(value);
      if (strcmp(param,"#zero") == 0)   ptr->zero = atoi(value);
  }
  fclose(fh_load);
  printf("Read numerical paramters:\nnx\t%4d\nny\t%4d\n", ptr->nx, ptr->ny);
  printf("lx\t%.4f\nly\t%.4f\n", ptr->lx, ptr->ly);
  printf("cfl number\t%.2f\nfinal time\t%.2f\n",ptr->cfl, ptr->t_max);
  printf("Starting from file\t%s\n", ptr->ic_name);
  //printf("Physical Parameters:\n");
  //printf("viscosity\t%.4e\n", ptr->visc);

  printf("Initial Conditions Parameters:\n");
  printf("U\t%.3f\na\t%.3f\n", ptr->U, ptr->a);
  printf("Artificial Viscosity Power = %d\n", 2*(ptr->p));

 
  printf("Forcing Parameters:\n");
  printf("forcing multiplier\t%f\n", ptr->FFlag);
  printf("Amplitude of Forcing %f\n", ptr->a_force);
  printf("Wavenumber of Forcing %f\n", ptr->k_force);

  ptr->vrt = malloc(sizeof(fftw_complex)*(ptr->nx)*(ptr->ny));
  ptr->vel.u = malloc(sizeof(fftw_complex)*(ptr->nx)*(ptr->ny));
  ptr->vel.v = malloc(sizeof(fftw_complex)*(ptr->nx)*(ptr->ny));
  ptr->grad_v.u = malloc(sizeof(fftw_complex)*(ptr->nx)*(ptr->ny));
  ptr->grad_v.v = malloc(sizeof(fftw_complex)*(ptr->nx)*(ptr->ny));
  ptr->sp = malloc(sizeof(fftw_complex)*(ptr->nx)*(ptr->ny));
  dx = (ptr->lx)/(ptr->nx); dy = (ptr->ly)/(ptr->ny);
  x = malloc(sizeof(fftw_complex)*(ptr->nx));
  y = malloc(sizeof(fftw_complex)*(ptr->ny));
  ptr->visc = pow((ptr->lx)/((ptr->alpha)*pi*(ptr->nx)),2*(ptr->p));
  printf("Set Viscosity Coefficient to %e\n", ptr->visc);
}
void load_text(ctrl_ptr ptr)
{
  char valuex[160];
  int j, k, s;
  j = 0; k = 0; s = 0;
  fh_load = fopen(ptr->ic_name, "r");
  if (fh_load == NULL) {
    printf("ic file missing\n");
    exit(1);
  }
  while (fgets(line, 160, fh_load) != NULL)
  {
      sscanf(line, "%s\n", valuex); 
      //printf("%s\t%f\n", valuex,atof(valuex));
      if (k < (ptr->ny)) 
      {
        ptr->vrt[k+(ptr->ny)*j] = atof(valuex);
        //printf("%d\t%d\t%e\n",j,k,ptr->vrt[k+(ptr->ny)*j]);
        k++;
      } else {
        k = 0; j++;
        ptr->vrt[k+(ptr->ny)*j] = atof(valuex);       
	//printf("%d\t%d\t%e\n",j,k,ptr->vrt[k+(ptr->ny)*j]);
	k++;
      }
      s++;
      if (s > (ptr->nx)*(ptr->ny)) {
        printf("Gridpoint number mismatch.\n");
        printf("Check that .in and ic_file have same number of gridpoints.\n");
        exit(1);
      }
  }
  fclose(fh_load);
}
void save_text(ctrl_ptr ptr, char *fname)
{
  int j, k;
  
  fh_save = fopen(fname,"w");
  fprintf(fh_save, "%%# 1.x 2.y 3.vorticity 4. u 5. v 6. grad_x vrt 7. grad_y vrt \n");
  fprintf(fh_save, "%%# 8. spectrum(vorticity) 9. kx 10. ky\n\n");
  find_spectrum(ptr);
  for (j = 0; j < (ptr->nx); j++){
    x[j] = -0.5*(ptr->lx)+(0.5+j)*dx;
    if (j > (ptr->nx)/2) kx = j - (ptr->nx);
    else kx = j;
    for (k = 0; k < (ptr->ny); k++) {
      if (k > (ptr->ny)/2) ky = k - (ptr->ny);
      else ky = k;
      y[k] = -0.5*(ptr->ly)+(0.5+k)*dy;
      fprintf(fh_save, "%e\t%e\t", x[j],y[k]);
      fprintf(fh_save, "%e\t%e\t", creal(ptr->vrt[k+(ptr->ny)*j]), creal(ptr->vel.u[k+(ptr->ny)*j]));
      fprintf(fh_save, "%e\t%e\t", creal(ptr->vel.v[k+(ptr->ny)*j]), creal(ptr->grad_v.u[k+(ptr->ny)*j]));
      fprintf(fh_save, "%e\t%e\t", creal(ptr->vel.u[k+(ptr->ny)*j]), cabs(ptr->sp[k+(ptr->ny)*j]) );
      fprintf(fh_save, "%e\t%e\n", kx, ky);
    }
    fprintf(fh_save,"\n");
  }
  fclose(fh_save);
}


